library(Seurat)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(Matrix)
library(cowplot)
## Warning: package 'cowplot' was built under R version 3.5.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
##
## ggsave
library(scrattch.hicat)
library(patchwork)
## Warning: package 'patchwork' was built under R version 3.5.3
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:cowplot':
##
## align_plots
#setwd('~/postdoc2/Gedankenpapers/CNevomanuscript/code/snRNAseq/mouse/processing/')
Load data.
load('../data/rawdata.RData')
p1<-DimPlot(object = N31, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = N31, reduction = 'tsne', label=TRUE,group.by = 'DCN')
p3<-DimPlot(object = N31, reduction = 'tsne', label=TRUE,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
plot number of detected genes and reads per cluster
VlnPlot(object = N31, feature = c("nFeature_RNA", "nCount_RNA"))
Plot maker genes
FeaturePlot(object = N31, features = c("Snap25","Slc17a6","Gad1","Slc6a5","Gabra6"), reduction= "tsne")
Subset data: clusters 1,0,18 gabra6+ clusters 4,9,14,15,17 low genes/counts
dcn<-subset(N31,idents=c(0,1,18,4,9,14,15,17),invert=T)
DimPlot(dcn)
dcn <- FindVariableFeatures(object = dcn,selection.method = "vst",nfeatures = 2000,verbose=FALSE)
dcn <- ScaleData(object = dcn,vars.to.regress = c('nCount_RNA','FACS','nFeature_RNA'))
## Regressing out nCount_RNA, FACS, nFeature_RNA
## Centering and scaling data matrix
#dcn <- RunPCA(object = dcn,npcs = 30,verbose = FALSE)
source('../../helperfunctions/pc_modification_functions_S3_forRNA.R')
## Loading required package: qlcMatrix
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'qlcMatrix'
dcn<-RunTruncatedPCA(dcn,n.genes.pc = 40)
## PC_ 1
## Positive: Zfhx3, Pdzd2, Grm7, Slc24a4, Fbxl7, Prkg1, Ptprk, Fxyd6, Foxp2, Ramp3
## Grm3, Asic2, Syt6, Coch, Trpm3, Gabrb1, Kcnip1, Fbln5, Cdh4, Gm40999
## Drd2, Grid2, Prr16, Ncald, Adarb2, Prex2, Glra3, Sorcs3, Trp63, Mal2
## Negative: Pde4b, Nfib, St18, Cdh20, Zfpm2, Gjc3, Zeb2, Rgs6, 1700047M11Rik, Ptprz1
## Tmtc2, Kit, Plxnb3, Mgat4c, Neat1, Npas3, Sec14l5, Sox2ot, Ttyh2, Cmtm5
## Mag, Slc6a5, Igsf11, Synpr, Cldn11, Fa2h, Tspan15, Stxbp2, Aspa, Cdh19
## PC_ 2
## Positive: Rgs6, Zfpm2, Mgat4c, Kit, Ptprz1, Nfib, Slc6a5, Rbfox1, Nxph1, Tmtc2
## Synpr, Npas3, Dab1, Slc36a2, Penk, Trim67, Hs3st4, Shisa9, Kcnip4, Plcb1
## Pde4b, Gm3764, Hs6st3, Stxbp2, Map3k8, Afap1l2, Kcnn3, Apcdd1, Hs3st2, Adrb2
## Negative: Ptprk, Grm7, Gjc3, Zfhx3, 1700047M11Rik, Neat1, Plxnb3, Sec14l5, Mog, Prr5l
## Sox2ot, Slc24a4, Mag, St6galnac3, Cdh19, Pdzd2, Cldn11, St18, Aspa, Prkg1
## Grm3, Fbxl7, Tspan2, ENSMUSG00000027199, Igsf11, Cmtm5, Foxp2, Bcas1, Mobp, Fxyd6
## PC_ 3
## Positive: Sec14l5, 1700047M11Rik, Aspa, Mog, St18, ENSMUSG00000111219, Mag, Cldn11, Mobp, Prr5l
## E330020D12Rik, Tnfaip6, Lpar1, Rasgrp3, Fa2h, Gm44866, Tspan15, Gas7, Nipal4, Hapln2
## Gjc3, Tnni1, Plxnb3, Cldn14, A230001M10Rik, Insc, Anln, Opalin, Apod, Rnf220
## Negative: Stk32a, Slc39a12, Slc38a3, Cd38, Slc4a4, Aqp4, Sox6, Etnppl, Acss3, Cyp2j9
## Gja1, Ntsr2, Gpr37l1, Pdgfd, Bmpr1b, Gjb6, A330076C08Rik, Aldh1a1, Gm6145, Acsbg1
## S1pr1, Slc7a10, Ednrb, F3, Rorb, Slc1a2, Pax3, Egfr, Serpine2, Agt
## PC_ 4
## Positive: Spp1, Cpne4, Unc5d, Kcnip4, Gpc5, D630045J12Rik, Rbfox1, Zfp804b, Nrg1, Grik1
## Fam19a2, Itpr1, Nell2, Adamts6, Zfyve28, Dab1, Kirrel3, Col24a1, Pld1, Chrna4
## Pde4b, Plpp3, Gm10801, Ninl, Zpbp, Gm26917, Slitrk6, Kank1, Gpr83, Epb41l2
## Negative: Adarb2, Shisa9, Slc6a5, Kit, Plcb1, Nxph1, Synpr, Syndig1, Ptprk, Grm3
## Asic2, Ptprz1, Slc24a4, Igfbp2, Kcnn3, Tmtc2, Fam107b, Ust, Slc36a2, Zfhx3
## Grid2, Apcdd1, Map3k8, Fxyd6, Plppr4, Syt6, Penk, Zfpm2, Ramp3, Nfib
## PC_ 5
## Positive: Kirrel3, Gria1, Zfp804b, Gm38048, Mgat4c, Fam19a2, Car8, Trabd2b, Ryr1, Sst
## Calb1, Prkg1, Mctp1, Unc5d, Gabrb1, Grik1, Grm7, Casq2, Nckap5, Gm28928
## Iltifb, Gabrg3, Stac, Cacng3, Myo10, Nell2, Cck, Col24a1, Fbn2, Prex2
## Negative: Spp1, Adarb2, Fbxl7, Pdzd2, Slc24a4, Stxbp2, Ramp3, Grm3, Zfhx3, Syt6
## Gm3294, Coch, Zfyve28, Itpr1, Tox2, Pld1, Plpp3, Aqp4, Gpr83, Etnppl
## Acsbg1, Itpr2, Zpbp, Hs6st3, Stk32a, Slc39a12, Slc7a10, Cd38, Gjb6, Qk
ElbowPlot(object = dcn,ndims = 30)
usefuldims=1:26
dims.remove=c(10,12,13,18,19,21,23,24,25)
usefuldims=usefuldims[!usefuldims %in% dims.remove]
dcn<- FindNeighbors(dcn,dims=usefuldims)
## Computing nearest neighbor graph
## Computing SNN
dcn <- FindClusters(object = dcn, reduction = "pca", dims = usefuldims, resolution = 2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 4988
## Number of edges: 201872
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6724
## Number of communities: 19
## Elapsed time: 0 seconds
dcn <- RunTSNE(object = dcn, dims = usefuldims, perplexity=30, dim.embed = 2,seed.use = 2)
p1<-DimPlot(object = dcn, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = dcn, reduction = 'tsne', label=TRUE,group.by = 'DCN')
p3<-DimPlot(object = dcn, reduction = 'tsne', label=TRUE,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
VlnPlot(dcn,features=c("nFeature_RNA", "nCount_RNA"))
Figure out which dimensions are actually useful for clustering.
#FeaturePlot(dcn, features=paste0("PC_", 1:30),cols = c("navy", "olivedrab1"))
plot maker genes
FeaturePlot(object = dcn, features = c("Snap25","Slc17a6","Gad1","Slc6a5","Gabra6"), reduction= "tsne")
let’s look only at the excitatory neurons.
ex<-subset(dcn,idents=c(13,6,10,4,5,8,1),invert=F)
DimPlot(ex)
ex <- FindVariableFeatures(object = ex,selection.method = "vst",nfeatures = 2000,verbose=FALSE)
ex <- ScaleData(object = ex,vars.to.regress = c('nCount_RNA','FACS','nFeature_RNA'))
## Regressing out nCount_RNA, FACS, nFeature_RNA
## Centering and scaling data matrix
ext<-RunTruncatedPCA(ex,dim=1:30)
## PC_ 1
## Positive: Grm7, Celf4, Prkg1, Kirrel3, Fnbp1l, Gsg1l, Dscaml1, March1, Snca, Kcnb2
## Syt16, Luzp2, Kcnq5, Chl1, Mgat4c, Slc6a1, Zfp804b, 6330403K07Rik, Cacna1e, Shisa9
## Unc5c, Sema6d, Atp2b4, Lmx1a, Fam19a2, Fxyd6, Dcc, Mal2, Pkib, Ablim3
## Negative: Spp1, Crhr1, Peak1, Kcnip4, Hs6st3, Stxbp2, Fign, Zfpm2, Galnt18, Fam110b
## Gpc6, Zfyve28, Adarb2, Tcerg1l, Pld1, Fignl2, L3mbtl4, Sgpp2, Gm26917, Zpbp
## Sgcd, Cpne4, Exog, Cnksr2, Slmap, Slc7a2, Opa1, Trpm2, Bmper, Anxa3
## PC_ 2
## Positive: Sorcs3, Cntn4, Kynu, Cacna1g, Fstl5, Sulf1, Cpne4, Ankfn1, March1, Gabrg3
## Fign, Necab1, Pde1a, Lamb1, Adamtsl1, Sst, Adcy2, Crhr1, Gal, Gm28928
## 6330403K07Rik, Galnt18, Fxyd6, Peak1, Sntg2, Tcerg1l, Prr16, Cpne6, Cacng3, Zfp521
## Negative: Cntnap4, Lmx1a, Kcnc2, Nrg1, Gm45904, Grm7, Thsd7b, Mgat4c, Fam19a2, Gpc5
## Zfp804b, Unc5c, Nox4, Nrxn3, Col24a1, Ndst4, Kcnip4, Prkg1, Alk, Gm11309
## Htr2a, Dab1, Gm13974, Kirrel3, Ecel1, Cd36, Penk, Grm3, Adamts6, Sema3e
## PC_ 3
## Positive: Nfib, Nxph1, Asic2, Prr16, Sgcd, Chrm2, Hs6st3, Gm32647, Gfra1, Dlgap2
## L3mbtl4, Glis3, Ptprk, Gm10800, Ano1, Ndst4, Ankfn1, Thsd7b, Gm10801, Penk
## Foxp2, Zeb2, Htr2a, Ust, Zfp521, Zfpm2, Sox6, Col24a1, Hs3st4, Syndig1
## Negative: Trpm3, Galntl6, Zfp804b, Atp1a3, Dscaml1, Cpne4, CT010467.1, Chl1, Scg2, Cpne6
## Syngr3, Aldoc, Cst3, Txndc15, Tcerg1l, Sntg2, Nox4, Fnbp1l, Mfsd5, Gsg1l
## Gpc6, Fign, Pld3, Ecel1, Ablim3, Tram1l1, Serpina9, Napa, P4hb, Gstm5
## PC_ 4
## Positive: Grik1, Trpm3, Cpne4, Ano3, Sulf1, Gm10801, Greb1, Gm10800, Lama2, Zfpm2
## Kcnc2, Nox4, Dab1, Cacna1e, Fam107b, Zfp804b, Gpc5, Chl1, Sema3e, Pde1a
## Syndig1, Trabd2b, Gm13974, Gm28928, Adcy2, Sdk1, Mgat4c, Prlr, Nckap5, Sox6
## Negative: Asic2, Atp1a3, CT010467.1, Lmx1a, Psap, Col24a1, Fxyd6, Gm32647, Thsd7b, Scg2
## Unc5d, Pkib, L3mbtl4, Nfib, Spp1, Atp2b4, Pld3, 6330403K07Rik, Crhbp, Mfsd5
## Aldoc, Cdh4, Serpina9, Mboat7, Luzp2, Tram1l1, Ndufs7, Napa, Cacybp, Rcn1
## PC_ 5
## Positive: Gm10800, Hpse2, Gm10801, Lmx1a, Lama2, Gm20754, Synpr, Prkg1, Scn3b, Gm32647
## Myo1b, Fxyd6, Atp2b4, Casd1, Bnc2, Slc18a2, Fbn2, Col24a1, Insc, Otoa
## Tpd52l1, Gm19965, Crhbp, Ush2a, Prrg1, Kcnn3, Gm42397, Asic2, Fbln1, Col5a2
## Negative: CT010467.1, Atp1a3, Grik1, Ano3, Nxph1, Trpm3, Dpp10, Kcnip4, Dcc, Afdn
## Prr16, Sez6l, Mctp1, Ndst4, Dab1, Sgpp2, Gfra1, Kcnb2, Sulf1, Ncoa2
## Cacna1e, Slc6a1, Cntn4, Sema6d, Adcy2, Sntg2, Nfib, Ankfn1, Zfp521, Sdk1
ElbowPlot(object = ext,ndims = 30)
#ext <- JackStraw(object = ext, dims=30)
#ext<- ScoreJackStraw(ext,dims=1:30)
#JackStrawPlot(object = ext,dims=1:30)
usefuldims=1:18
dims.remove=c(17)
dims.remove=c()
usefuldims=usefuldims[!usefuldims %in% dims.remove]
replot with dropped pcs
ext<- FindNeighbors(ext,dims=usefuldims)
## Computing nearest neighbor graph
## Computing SNN
ext <- RunTSNE(object = ext, dims = usefuldims, perplexity=30, dim.embed = 2)
ext <- FindClusters(object = ext, reduction.type = "pca", dims.use = usefuldims, resolution = 2, print.output = 0, save.SNN = TRUE, force.recalc = TRUE)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2298
## Number of edges: 88158
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6727
## Number of communities: 18
## Elapsed time: 0 seconds
p1<-DimPlot(object = ext, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = ext, reduction = 'tsne', label=TRUE,group.by = 'DCN')
p3<-DimPlot(object = ext, reduction = 'tsne', label=TRUE,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
merge clusters using allen package.
rd.dat <- t(GetAssayData(object = ext, slot = "scale.data"))
merge.param <- de_param(padj.th = 0.05,
lfc.th = 1,
low.th = 1,
q1.th = 0.4,
q.diff.th = 0.5,
de.score.th = 40)
merge.result <- merge_cl(as.matrix(GetAssayData(object = ext, slot = "data")),
cl = ext$RNA_snn_res.2,
rd.dat = rd.dat,
de.param = merge.param)
## Loading required package: limma
ext$merged.res.2<-as.factor(paste('ex_',merge.result$cl,sep=''))
p2<-DimPlot(ext,group.by = "merged.res.2",label=T)
p1<-DimPlot(ext,group.by = "RNA_snn_res.2")
p3<-DimPlot(ext,group.by = 'DCN')
plot_grid(p1,p2,p3,ncol=3)
Idents(ext)<-'merged.res.2'
p2<-DimPlot(ext,group.by = "merged.res.2")
p1<-DimPlot(ext,group.by = "FACS")
p3<-DimPlot(ext,group.by = 'DCN')
plot_grid(p2,p3,p1,ncol=3)
#save(file='ex_celltypes.RData',ext)
now look at non-glycinergic inhibitory cell types
inh<-subset(dcn,idents=c(0,2,9),invert=F)
DimPlot(inh)
inh <- FindVariableFeatures(object = inh,selection.method = "vst",nfeatures = 2000,verbose=FALSE)
inh <- ScaleData(object = inh,vars.to.regress = c('nCount_RNA','FACS','nFeature_RNA'))
## Regressing out nCount_RNA, FACS, nFeature_RNA
## Centering and scaling data matrix
source('../../helperfunctions/pc_modification_functions_S3_forRNA.R')
## Loading required package: qlcMatrix
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'qlcMatrix'
inh<-RunTruncatedPCA(inh,n.genes.pc = 40)
## PC_ 1
## Positive: Ptprt, Malat1, Lgr5, Angpt1, Cd226, Arhgap15, Mgat4c, Cryzl2, Npsr1, Kcnh8
## ENSMUSG00000039840, Abi3bp, Pik3c2g, Slco4a1, Arhgap28, Cdc14a, Esr1, Pde4b, Spata13, Nectin3
## Adcy2, Ust, Gm28928, Cep128, Cpne8, Btk, Prdm1, Chst7, Gm15983, Sncaip
## Negative: CT010467.1, Ndrg4, Syvn1, Gm26917, Adcy1, Peg3, Psap, Mdm4, Tm9sf2, Map2k5
## Arsa, Ncald, Appl2, Snph, Sez6l2, Pisd, Entpd6, Pfkfb3, Dhcr7, Slc25a1
## Clcn7, Cyp51, Hsd17b7, Slc4a3, Hivep1, Timm50, Rac1, Lnx1, Gpaa1, Naxd
## PC_ 2
## Positive: Tenm2, Grm8, CT010467.1, Cadps2, Cacng3, Igfbp5, Gfra1, Anxa11, Bscl2, Bves
## Igfbp4, Mgat4c, Rac1, Pigt, Stc1, Slc19a1, Slc15a2, Odf3b, Pmm1, Pde4b
## Adgrg1, Serpinf1, Slc25a39, Disp3, Arhgef26, Arhgap28, Cptp, Pcsk2, Thumpd1, Mpp6
## Negative: Kcnc2, Plcl1, Klhl29, Rbfox1, Kcnip4, Ccbe1, Prkca, Shisa9, Galnt18, Gm35721
## Fam19a2, Plekha7, Prkd1, Myo18b, Rxfp1, Slit3, Grin2a, Neurl1b, Map2k5, Cyp4f18
## Rspo3, Pmfbp1, Fam81a, Lrmda, Lrriq1, Il1rapl2, Parn, Ly96, 2310002F09Rik, Pfkfb3
## PC_ 3
## Positive: Ntng1, Rmst, Gfra1, Tenm2, Taf13, Tor3a, Rbfox1, Ptprt, Fam19a2, Nubpl
## Mgat4c, Arhgef26, Il1rapl2, Grm8, Chn2, Grin2a, Adcy2, Irf5, Vcl, Gramd1c
## Wipf1, Gm15601, D630003M21Rik, Trp53rka, Dlgap2, Lrmda, Smc2, Pde4b, Pld5, Gm15983
## Negative: Klhl29, Lax1, Chrna4, Slc7a3, Tor1a, Fam163b, Kcnc2, Gstm1, Rpl22, Fam192a
## Slc35b1, Mfsd3, Prkd1, Tdrkh, Gnb3, Doc2a, D930016D06Rik, Nkain4, Ada, Ppp4r1l-ps
## Asrgl1, Ubl3, Lss, Pgd, Xkr8, Ninj1, Znhit6, Btbd19, Golga7b, Kcnip4
## PC_ 4
## Positive: Il1rapl2, Rmst, Fam19a2, Tenm2, Kcnip4, Kcnc2, Serpinf1, Cadps2, Rbfox1, Rxfp1
## Ctnna1, Syvn1, Myo18b, Cdh11, Nprl3, Gstm4, 4931423N10Rik, Fam163b, Gm12371, Col4a2
## Mcfd2, Mgat4c, Sf3b3, Fam192a, Ifngr1, Stat2, Ntng1, 1810058I24Rik, Tonsl, 4931406C07Rik
## Negative: Gm32647, Klhl29, 2900079G21Rik, Myt1, Tmem161a, Rcc1, Arhgef26, Tefm, Tbc1d4, Tars
## Spata9, Anxa11, Zmiz1os1, Hmgxb3, Pantr1, Gprin1, Cngb1, Dlgap2, Naf1, Mfsd5
## Doc2a, 1700102P08Rik, Pik3r6, Caskin1, Myh9, Snph, Slc5a5, Drd2, Smpdl3a, Dcbld2
## PC_ 5
## Positive: C230088H06Rik, Gfra1, Tfdp1, 2310002F09Rik, Galnt18, Malat1, D030028A08Rik, Gm32647, Pmfbp1, Arhgef26
## Zfp872, Klhl29, Lnx1, Hic2, 4931415C17Rik, Timmdc1, Ddx3x, Mad2l1, Gm43848, Xrcc5
## Mpp6, Otulin, Uap1l1, Fam210a, Camkk1, Airn, Spry4, Fto, Cxxc1, Adat1
## Negative: Ly96, Gm17949, Nek3, Mfsd12, Serpinf1, Ccdc65, Tprkb, Ret, Opn3, Smpdl3a
## Serpine1, CT010467.1, Tm4sf1, ENSMUSG00000015980, Btk, Dlec1, Tubg1, Gm35721, Il21r, Msto1
## Cryzl2, Mettl1, Gm6859, Gm830, D530049I02Rik, Adcy2, Chrna4, Smyd5, Rdh11, Spp1
ElbowPlot(object = inh)
#inh <- JackStraw(object = inh, dims=20)
#inh <- ScoreJackStraw(inh,dims=1:20)
#JackStrawPlot(object = inh,dims=1:20)
usefuldims=1:4
dims.remove=c(1,3)
usefuldims=usefuldims[!usefuldims %in% dims.remove]
inh<- FindNeighbors(object = inh,dims=usefuldims)
## Computing nearest neighbor graph
## Computing SNN
inh<- FindClusters(object = inh, resolution = 2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1401
## Number of edges: 31880
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7305
## Number of communities: 22
## Elapsed time: 0 seconds
inh <- RunTSNE(object = inh, dims = usefuldims, perplexity=30, dim.embed = 2)
p1<-DimPlot(object = inh, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = inh, reduction = 'tsne', label=TRUE,group.by = 'DCN')
p3<-DimPlot(object = inh, reduction = 'tsne', label=TRUE,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
#FeaturePlot(inh, features=paste0("PC_", 1:4),cols = c("navy", "olivedrab1"))
merge clusters using allen package.
rd.dat <- t(GetAssayData(object = inh, slot = "scale.data"))
merge.param <- de_param(padj.th = 0.05,
lfc.th = 1,
low.th = 1,
q1.th = 0.4,
q.diff.th = 0.6,
de.score.th = 40)
merge.result <- merge_cl(as.matrix(GetAssayData(object = inh, slot = "data")),
cl = inh$RNA_snn_res.2,
rd.dat = rd.dat,
de.param = merge.param)
if (is.null(merge.result))
{inh$merged.res.2<-1
} else {
inh$merged.res.2<-merge.result$cl
}
p2<-DimPlot(inh,group.by = "merged.res.2")
p1<-DimPlot(inh,group.by = "RNA_snn_res.2")
p3<-DimPlot(inh,group.by = 'DCN')
p4<-DimPlot(inh,group.by = 'FACS')
plot_grid(p2,p3,p4,ncol=3)
Idents(inh)<-"merged.res.2"
VlnPlot(inh,features=c("nFeature_RNA", "nCount_RNA"))
plot marker genes
FeaturePlot(object = inh, features = c("Snap25","Slc17a6","Gad1","Slc6a5"), reduction= "tsne")
Now, lets do this for the big group of glycinergic cells
gly<-subset(dcn,idents=c(3,7,15),invert=F)
DimPlot(gly)
gly <- FindVariableFeatures(object = gly,selection.method = "vst",nfeatures = 2000,verbose=FALSE)
gly <- ScaleData(object = gly,vars.to.regress = c('FACS','nFeature_RNA'))
## Regressing out FACS, nFeature_RNA
## Centering and scaling data matrix
source('../../helperfunctions/pc_modification_functions_S3_forRNA.R')
## Loading required package: qlcMatrix
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'qlcMatrix'
gly<-RunTruncatedPCA(gly,n.genes.pc = 40)
## PC_ 1
## Positive: Cadm1, Adamtsl1, Spon1, Kirrel3, Trpm3, Ly6d, Syt16, Slc24a3, Unc5d, Cacna1g
## Fbn2, Megf10, Dab1, Shisa6, Fxyd6, Cdh1, Gabbr2, Igsf8, Rftn1, Tunar
## Gpc5, mt-Rnr2, Eps8, Arhgef28, Plekhg1, Sema4a, Plch2, Cnpy1, Nkain1, Ptpro
## Negative: Slc6a5, Kcnip4, Adam23, Penk, Ccbe1, Ano4, Rbms1, Clstn2, Tenm4, Mdh1
## Ndrg4, Slc6a12, Adamts3, ENSMUSG00000013539, Gm26917, Ttll5, Lpcat2, Fggy, Ptar1, Ablim2
## Rab31, Srgap2, Hdlbp, Nova1, Rbfox1, Arhgef6, Epb41l4b, Fhit, Mrps9, Lrmda
## PC_ 2
## Positive: Slc6a5, Mdh1, Kcnip4, Ccbe1, Itpr3, Map3k21, Penk, Stmn3, Cplx2, Mettl7a1
## Tgfb2, Ttr, Gjb1, Gm26597, Tbl3, Gm25127, 2610307P16Rik, Ctns, Rpusd4, Gm20609
## Ctxn1, Fyb, Trim28, Tmem53, Bysl, AC156026.1, Gm38155, Kcng2, Gm12240, Gm16283
## Negative: Cadm1, Spon1, Kirrel3, Slc24a3, Dab1, Plekhg1, Tnr, Ly6d, Gabbr2, Osbpl10
## Cacna1g, Adamts17, Atxn2l, Shisa6, Megf10, Ptprn2, Mical3, Plxna4, Anapc5, Sgcd
## Ptch1, Unc5d, Ptpro, Tcerg1l, Adamtsl1, Leng8, Trpm3, Syt16, Nkain1, Trim67
## PC_ 3
## Positive: CT010467.1, Atp1a3, Ndrg4, Stmn3, Ppp2r1a, B4gat1, Psap, Syngr3, Kcnn3, Mdh1
## Ckmt1, Cnbp, Dab1, Chpf, Ddost, Rrp1, Sqstm1, Pcbp1, Disp2, Slc39a6
## Syvn1, Reep5, Paf1, Atxn2l, Nlgn2, Yif1a, Tmem109, Dcaf7, Adora1, Tunar
## Negative: Tenm4, Slit2, Cntnap5c, Asic2, Kcnip4, Bicc1, Rorb, Zfhx3, Zeb2, Cntn5
## 2610307P16Rik, Il1rapl2, Grin2a, Rab27b, Gpc5, Sgcd, Rsu1, Chst8, Lrmda, Diaph2
## Gm2164, Gm20449, Ctnna3, Map3k21, Gpr161, Slc10a7, Rbfox1, Lpcat2, Wwc1, Efl1
## PC_ 4
## Positive: Npas3, Cntnap5a, Dab1, Pkd2, C2cd3, Afap1, Ptprn2, Osbpl3, Etv5, Eml4
## Rbfox1, Raph1, Pip4k2a, Fam193a, Manba, Tctn1, Sema4a, Limd2, Wscd1, ENSMUSG00000022253
## Asap2, Rnf146, ENSMUSG00000000325, Inpp4a, Nkain1, Fhit, Trpc7, Yaf2, Ptp4a2, Chd7
## Negative: Tenm3, Cntnap5c, Pdzrn4, Chsy3, Ebf3, Meis2, Fam19a1, Cntn5, Pdgfb, CT010467.1
## Gm6999, Tmem132e, Trim36, Mfsd5, Etnppl, Tle2, Spink4, Clec2j, Card10, Tst
## 4921539H07Rik, Gm15518, Mdh1, C430002N11Rik, C230012O17Rik, Atp1a3, Syngr3, Mir1969, Ebf1, Asic2
## PC_ 5
## Positive: Asic2, Arhgap6, Zeb2, Pde8b, Tenm4, Dnhd1, Ccbe1, Grin2a, Ssbp3, Slco4a1
## Cdh6, Rab3b, Plch2, Tmem28, Tmem74b, Lin9, Xpo5, Tmem231, Syt16, Rorb
## Fam131a, Phf12, Errfi1, Pnpla2, Cabyr, Slc6a5, Capns1, Siah3, Paf1, Zfhx3
## Negative: Cntnap5a, Tst, Cntnap5c, Il1rapl2, Dab1, Chsy3, Ebf3, Spp1, Kcnn3, Etnppl
## Rnf17, Frmd8os, Thsd7b, Fastkd2, Gm6999, Eva1a, Mir1969, Srgap1, Dcst1, Rbfox1
## Gm15518, 4921539H07Rik, Tph1, Tle2, C230012O17Rik, Fam160a2, Ankzf1, Glb1, Gm14696, Rny1
ElbowPlot(object = gly)
gly <- JackStraw(object = gly, dims=20)
gly <- ScoreJackStraw(gly,dims=1:20)
JackStrawPlot(object = gly,dims=1:20)
usefuldims=1:5
dims.remove=c()
usefuldims=usefuldims[!usefuldims %in% dims.remove]
gly<- FindNeighbors(object = gly,dims=usefuldims)
## Computing nearest neighbor graph
## Computing SNN
gly<- FindClusters(object = gly, reduction.type = "pca", dims.use = usefuldims, resolution = 2, print.output = 0, save.SNN = TRUE, force.recalc = TRUE,algorithm = 1,n.start = 20)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 851
## Number of edges: 26235
##
## Running Louvain algorithm...
## Maximum modularity in 20 random starts: 0.5111
## Number of communities: 13
## Elapsed time: 0 seconds
gly <- RunTSNE(object = gly, dims = usefuldims, perplexity=30, dim.embed = 2)
p1<-DimPlot(object = gly, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = gly, reduction = 'tsne', label=TRUE,group.by = 'DCN')
p3<-DimPlot(object = gly, reduction = 'tsne', label=TRUE,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
merge clusters using allen package.
rd.dat <- t(GetAssayData(object = gly, slot = "scale.data"))
merge.param <- de_param(padj.th = 0.05,
lfc.th = 1,
low.th = 1,
q1.th = 0.4,
q.diff.th = 0.6,
de.score.th = 40)
merge.result <- merge_cl(as.matrix(GetAssayData(object = gly, slot = "data")),
cl = gly$RNA_snn_res.2,
rd.dat = rd.dat,
de.param = merge.param)
if (is.null(merge.result))
{gly$merged.res.2<-1
} else {
gly$merged.res.2<-merge.result$cl
}
p2<-DimPlot(gly,group.by = "merged.res.2")
p1<-DimPlot(gly,group.by = "RNA_snn_res.2")
p3<-DimPlot(gly,group.by = 'DCN')
p4<-DimPlot(gly,group.by = 'FACS')
plot_grid(p2,p3,p4,ncol=3)
#FeaturePlot(gly, features=paste0("PC_", 1:10),cols = c("navy", "olivedrab1"))
Idents(gly)<-"merged.res.2"
Now let’s integrate these blobs of inhibitory or glycinergic cells with the smaller distinct groups of those cells and save out the cell types.
gly.complete<-subset(dcn,idents = c(3,7,15,12,16))
#annotate
prelim.clusters<-as.data.frame(gly.complete$RNA_snn_res.2)
prelim.clusters$`gly.complete$RNA_snn_res.2`<-as.character(prelim.clusters$`gly.complete$RNA_snn_res.2`)
prelim.clusters[names(gly$merged.res.2),]<-paste('gly_',as.character(gly$merged.res.2))
prelim.clusters$`gly.complete$RNA_snn_res.2`<-as.factor(prelim.clusters$`gly.complete$RNA_snn_res.2`)
colnames(prelim.clusters)<-'prelim.clusters'
gly.complete$prelim.clusters<-(prelim.clusters)
Idents(gly.complete)<-'prelim.clusters'
DimPlot(gly.complete)
gly.complete <- FindVariableFeatures(object = gly.complete,selection.method = "vst",nfeatures = 2000,verbose=FALSE)
gly.complete <- ScaleData(object = gly.complete,vars.to.regress = c('nCount_RNA','FACS','nFeature_RNA'))
## Regressing out nCount_RNA, FACS, nFeature_RNA
## Centering and scaling data matrix
gly.complete<-RunTruncatedPCA(gly.complete,n.genes.pc = 40)
## PC_ 1
## Positive: Penk, Ccbe1, Jup, Slc36a2, Npas3, Trim67, Stxbp2, Etv5, Shisa9, Apcdd1
## Cdh6, Runx2, Lrrn1, Gfra2, Lpcat2, Slc6a12, Slc6a5, Ptar1, Atp1a3, Ndrg4
## Adrb2, Aqp6, Osbpl3, Tbc1d9, Hs3st4, Zyx, Tmem132a, Stmn3, Nuak2, Hs3st2
## Negative: Grm7, Ptprk, Unc5d, Edil3, Fstl5, Egfem1, Gm26871, Tenm3, Dpp10, Gria1
## Alcam, Dcc, Alk, Ebf3, Pcdh9, Gabrg3, Cacng3, Grik1, Gabrb1, Nell1
## Cntn5, Fam19a2, Kcnj6, Cadm1, Nckap5, Grm8, Megf11, Kcnk9, Sst, Frmpd3
## PC_ 2
## Positive: Shisa9, Nckap5, Gabrg3, Megf11, Gria1, Gm20754, Sst, Alk, Slc16a2, Cacng3
## Edil3, Kcnh7, Arhgap6, Cnih3, Cdh4, Plpp4, Slc30a3, Alcam, Syt16, Slc6a11
## Gap43, Pxdn, Fam126a, Slit1, Fam107b, Dcc, Gabrb1, Litaf, Gpd1l, Cidea
## Negative: Grik1, Slit2, Fam19a2, Zfhx4, Pde9a, Cntnap5c, AC124532.2, Grm8, Kcnk9, Egfem1
## Lrfn2, Tenm3, Adamts17, Nell1, Cacna1i, Grm4, Il1rapl2, Olfml2b, B3gat1, Zfp804b
## Skap2, Colq, Kcng4, Hydin, Lrmda, St6gal1, Sema3e, D630045J12Rik, Gpc5, Hhatl
## PC_ 3
## Positive: Pcdh15, Slc6a5, Tenm4, Kcnip4, Syt1, Alk, Ptprk, Megf11, Adam23, Gabrg3
## Edil3, Pcdh9, Rbfox1, Ror1, Gm20754, Cntn5, Asic2, Dcc, Adamts3, Gpc6
## Dlg2, Ttll5, Gria1, Kcnh7, Ebf3, Penk, Slit3, Cacng3, Ccbe1, Ppargc1b
## Negative: Dab1, Spon1, Cadm1, Ly6d, Kirrel3, Cdh1, Tspan18, Calca, Tox, Arhgef28
## Syt16, Adamtsl1, Megf10, Eps8, Nkain1, mt-Rnr2, Chrna6, Pde1a, Slc6a2, Th
## Syt5, Fbn2, Cacna2d1, Tgfa, Slc7a3, Wdr72, Wls, Fxyd6, Ptpro, Slc16a12
## PC_ 4
## Positive: CT010467.1, Atp1a3, Ndrg4, Ppp2r1a, Syngr3, Atxn2l, Ckmt1, Leng8, Hnrnpk, B4gat1
## Fxyd6, Disp2, Slc30a3, Dab1, Reep5, Acp2, Psap, Cnbp, Sgta, Fam163b
## Cars, Kcnn3, Stmn3, Trim67, Dcaf7, Nlgn2, Rab34, Anapc5, Syp, Tmub2
## Negative: Map3k21, Wdr72, Pcdh15, Gpc6, 2610307P16Rik, Litaf, Zfhx3, Afap1, Kcnip4, Zfp804b
## Rxfp1, Fgd5, Gm15477, Prr5l, Syt1, Npas3, Dlg2, Slc38a3, Grm7, Tenm4
## Dcc, Mettl7a1, Hydin, 4430402I18Rik, Alk, Slit1, Scn7a, ENSMUSG00000032648, 2310039L15Rik, Lamc2
## PC_ 5
## Positive: Tspan18, Wls, Th, Slc6a2, Chrna6, Cdh8, Tox, Calca, Cacna2d1, Dubr
## Disp3, Ctxn1, Fam159b, Irgm1, Kcnip4, Fam163b, Fibcd1, Plin5, Gpx3, Vwa5b1
## Lamc2, Chrna3, 4931415C17Rik, Cntnap5c, Gm45592, Zc2hc1c, Pcdh15, Gm40348, Wdr72, Gm5455
## Negative: Kirrel3, Ly6d, Cadm1, Adamtsl1, Dab1, Trpm3, Spon1, Megf10, Wscd1, Cnpy1
## Nol4, Pcdh9, Ptprk, Cdh1, Plekhg1, Inpp5a, Dpp10, Adamts17, Gpc5, Ebf3
## Nkain1, Sst, Tcerg1l, Fbn2, Plpp4, Gabrg3, Sema4a, Sergef, Shc4, Rin2
ElbowPlot(object = gly.complete)
usefuldims=1:10
dims.remove=c()
usefuldims=usefuldims[!usefuldims %in% dims.remove]
gly.complete<- FindNeighbors(object = gly.complete,dims=usefuldims)
## Computing nearest neighbor graph
## Computing SNN
gly.complete <- RunTSNE(object = gly.complete, dims = usefuldims, perplexity=30, dim.embed = 2)
p1<-DimPlot(object = gly.complete, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = gly.complete, reduction = 'tsne', label=TRUE,group.by = 'DCN')
p3<-DimPlot(object = gly.complete, reduction = 'tsne', label=TRUE,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
rename idents
gly.complete <- RenameIdents(object = gly.complete, `gly_ 0` = "gly_0")
gly.complete <- RenameIdents(object = gly.complete, `gly_ 1` = "gly_1")
gly.complete <- RenameIdents(object = gly.complete, `gly_ 4` = "gly_2")
gly.complete <- RenameIdents(object = gly.complete, `12` = "gly_3")
gly.complete <- RenameIdents(object = gly.complete, `16` = "gly_4")
gly.complete[["prelim.clusters"]] <- Idents(object = gly.complete)
p1<-DimPlot(object = gly.complete, reduction = 'tsne', label=TRUE)
p2<-DimPlot(object = gly.complete, reduction = 'tsne', label=F,group.by = 'DCN')
p3<-DimPlot(object = gly.complete, reduction = 'tsne', label=F,group.by = 'FACS')
plot_grid(p1,p2,p3,ncol = 3)
FeaturePlot(gly.complete,c("Slc17a6","Gad1","Slc6a5","Slc6a1"))
Now let’s port the subclustering results back to the integrative dcn dataset
final.clusters<-as.data.frame(dcn$RNA_snn_res.2)
final.clusters$`dcn$RNA_snn_res.2`<-as.character(final.clusters$`dcn$RNA_snn_res.2`)
final.clusters[names(ext$merged.res.2),]<-paste('',as.character(ext$merged.res.2),sep='')
final.clusters[names(inh$merged.res.2),]<-paste('inh_',as.character(inh$merged.res.2),sep='')
final.clusters[names(gly.complete$prelim.clusters),]<-paste('',as.character(gly.complete$prelim.clusters),sep='')
final.clusters$`dcn$RNA_snn_res.2`<-as.factor(final.clusters$`dcn$RNA_snn_res.2`)
colnames(final.clusters)<-'final.clusters'
dcn$final.clusters<-(final.clusters)
p1<-DimPlot(dcn,group.by = "final.clusters",pt.size = 1,label = T)
p2<-DimPlot(dcn,group.by = 'DCN',pt.siz=1)
p3<-DimPlot(dcn,group.by = 'FACS')
plot_grid(p1,p2,ncol = 2)
FeaturePlot(dcn,c("Snap25","Slc17a6","Gad1","Slc6a5"))
Idents(dcn)<-'final.clusters'
VlnPlot(dcn,c('nFeature_RNA','nCount_RNA'),ncol = 1)